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Abstract 

Molecular phenotypes are important links between genomic information and or- 
ganismic functions, fitness, and evolution. Complex phenotypes, which are also called 
quantitative traits, often depend on multiple genomic loci. Their evolution builds on 
genome evolution in a complicated way, which involves selection, genetic drift, mu- 
tations and recombination. Here we develop a coarse-grained evolutionary statistics 
for phenotypes, which decouples from details of the underlying genotypes. We derive 
approximate evolution equations for the distribution of phenotype values within and 
across populations. This dynamics covers evolutionary processes at high and low re- 
combination rates, that is, it applies to sexual and asexual populations. In a fitness 
landscape with a single optimal phenotype value, the phenotypic diversity within pop- 
ulations and the divergence between populations reach evolutionary equilibria, which 
describe stabilizing selection. We compute the equilibrium distributions of both quan- 
tities analytically and we show that the ratio of mean divergence and diversity depends 
on the strength of selection in a universal way: it is largely independent of the pheno- 
type's genomic encoding and of the recombination rate. This establishes a new method 
for the inference of selection on molecular phenotypes beyond the genome level. We 
discuss the implications of our findings for the predictability of evolutionary processes. 

1 Introduction 

In recent years, we have witnessed an enormous growth of information from genome sequence 
data, which has enabled large-scale comparative studies within and across species. How this 
genomic information translates into biological functions is much less known. Molecular 
functions integrate the genomic information of their constitutive sites, and they can often 
be associated with specific phenotypes. Many such phenotypes are quantitative traits: they 



1 



have a continuous spectrum of values and depend on multiple genomic sites. For example, 
the binding of a transcription factor to a regulatory DNA site can be monitored by its effect 
on the expression level of the regulated gene. 

In an evolutionary context, biological functions and their associated pheno types are quan- 
tified by their contribution to the fitness of an organism. Such fitness effects can in principle 
be uncovered from genomic data by comparative analysis. However, a genome-based analysis 
of phenotypic evolution can be exceedingly complicated. The source of these complications 
is two-fold: Quantitative traits often depend on numerous and in part unknown genomic 
sites. Moreover, the evolution of these sites is coupled by fitness interactions (epistasis) and 
by genetic linkage. For both reasons, the genomic basis of a complex phenotype is not com- 
pletely measurable. At the same time, details of site content, linkage, and epistasis should 
not matter for the evolution of the phenotype itself. This calls for an effective, coarse-grained 
picture of the evolutionary process at the phenotypic level, which is the topic of this paper. 
We will show that complex quantitative traits have universal phenotypic observables, which 
decouple from the trait's genomic basis. Such universality turns out to be important for the 
practical analysis of a quantitative trait: it provides a way to infer its fitness effects based 
solely on phenotypic measurements. 

The map from genotype to phenotype is a challenging problem for statistical theory. 
The reason is that epistasis and linkage generate correlations in a population: the popu- 
lation frequency of individuals with a combination of alleles at a set of genomic sites may 
be larger or smaller than the product of the single-site allele frequencies. We refer to these 
correlations by the standard term linkage disequilibrium (which is quite misleading, be- 
cause linkage correlations have nothing to do with disequilibrium). Linkage disequilibrium 
is strongest in asexually evolving populations, but it can also be maintained under sexual 
reproduction, whenever recombination between genomic loci is too slow to randomize allele 
associations [12]. Linkage disequilibrium and epistasis can make the genomic evolution of a 
quantitative trait a strongly correlated many-"particle" process, and these correlations are 
crucial for the resulting phenotype statistics. Any quantitative understanding of this dy- 
namics must be based on an evolutionary model that contains selection, mutations, genetic 
drift, and (in sexual populations) recombination - and is yet analytically tractable at least 
in an approximate way. Before we turn to the agenda of this paper, we briefiy summarize 
current models of genome evolution and their application to quantitative traits. 

All known analytically solvable genome evolution models for multiple sites are based on 
the assumption that linkage correlations vanish [5,17,20,33-35,57] or are small [6,42]. There 
are two classes of quantitative traits to which these models can be applied. One of these 
consists of phenotypes that depend only on a small number of genomic sites. Such phenotypes 
are mostly monomorphic and occasionally polymorphic at a single of their constitutive sites. 
Hence, allele changes at different sites are well separated in time and linkage disequilibrium 
is small, regardless of the level of recombination. An example of such microscopic traits 
is transcription factor binding sites in prokaryotes and simple eukaryotes, which typically 
have about ten functional bases. In a time-independent fitness landscape, the genomic and 
phenotypic evolution of microscopic traits leads to simple equilibrium states of Boltzmann 
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form, which can be used for the inference of selection (this type of equihbrium is reviewed 
in the next section) [7,49]. In this way, fitness landscapes for transcription factor binding, 
which depend on the binding energy as molecular phenotype, have been inferred from site 
sequence data in bacteria and yeast [38,39]. 

The other, complementary class is phenotypes with a large number of constitutive sites 
which are assumed to evolve under rapid recombination, so that linkage correlations remain 
small. This assumption is justified in sexual populations, if all of the sites are at suffi- 
cient sequence distance from each other. It implies that the phenotype distribution in a 
population is completely determined by the allele frequencies at the constitutive sites. Ex- 
amples are an organism's height, complex disease phenotypes or longevity, which depend on 
multiple genes on different chromosomes. Such macroscopic traits are always polymorphic 
at multiple constitutive sites, which leads to a distribution of trait values in a population. 
Macroscopic traits are the traditional subject of quantitative genetics, which focuses on a 
phenomenological description of these trait distributions [1,5,9,20,22,27,34,35,46,56,57]. 
Rapid recombination is a crucial ingredient for existing evolutionary models of macroscopic 
traits [17,33]. As long as linkage disequihbrium remains small, macroscopic traits also reach 
genomic and phenotypic Boltzmann equilibria in a time-independent fitness landscape (for 
details, see next section) [16,17,22,33,55,57]. 

Many interesting molecular phenotypes, however, cannot be assumed to evolve close to 
hnkage equihbrium. The stability of protein and RNA folds depend on their coding se- 
quence [21,51], protein binding affinities depend on the nucleotides encoding the binding 
domain [7,8], complex regulatory interactions depend on cis- regulatory modules with several 
binding sites [15,44], histone-DNA binding involves segments of about 150 base pairs [45]: 
these are typical examples of intermediate-level phenotypes with tens to hundreds of constitu- 
tive DNA sites. Such mesoscopic phenotypes, which are often building blocks of macroscopic 
traits, are generically polymorphic at several constitutive sites. In asexual populations, meso- 
scopic traits always evolve under substantial linkage disequilibrium. This dynamics governs, 
for example, the evolution of antibiotic resistance in bacteria [53] and the antigenic evolution 
in human influenza A [52] . But mesoscopic traits can build up linkage disequilibrium even in 
sexual populations, because their constitutive sites are localized in a small genomic region, 
which limits the power of recombination [12]. Genomic evolution of multiple sites under 
weak recombination is a strongly correlated process, which generates cooperative phenom- 
ena such as clonal interference and background selection [2,11,18,24,43,47,48]. In other 
words, the phenotype distribution in a population is no longer determined by the allele fre- 
quencies at the constitutive sites, but it depends on the full distribution of genotypes. In 
this volume, Shraiman and colleagues show that the buildup of sequence correlations with 
decreasing recombination rate leads to a transition from allele selection to genotype selection, 
which is analogous to the glass transition in the thermodynamics of disordered systems [50]. 
These correlations lead to the breakdown of known analytical models for quantitative trait 
evolution. 

The evolution of molecular traits under genetic hnkage is the focus of this paper. Our 
dynamical model for phenotypes is grounded on the evolution of their constitutive genotypes 
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by selection, mutations, and genetic drift, which is reviewed in Section 2. In Section 3, we 
derive approximate, self-consistent equations for the asexual evolution of trait values in a 
population, which are parametrized by their mean and variance (called trait diversity). We 
show that this dynamics is quite different from trait evolution in sexual populations: In a 
time-independent fitness landscape, the joint distribution of trait mean and diversity con- 
verges to a non-equilibrium stationary state, yet the marginal distributions of both quantities 
still reach solvable evolutionary equilibria. In Section 4, we apply this model to evolution 
in a fitness landscape with a single trait optimum, where these equilibria describe the trait 
statistics under stabilizing selection. We compute the expected equilibrium diversity within 
populations, the divergence across populations, and the distance of a population from the 
fitness peak. In Section 5, we derive the statistics of population fitness and entropy in the 
equilibrium ensemble. Specifically, we compute the genetic load, which is defined as the 
difference between the maximum fitness and the mean population fitness, and the fitness 
fiux, which quantifies the total amount of adaptation between the neutral state and the 
stabilizing-selection equilibrium. The equilibrium entropy statistics is shown to determine 
the predictability of evolutionary processes from single-population data. All our analytical 
results are confirmed by numerical simulations. 

Throughout the paper, we compare our derivations and results for non-recombining pop- 
ulations with their counterparts for rapid recombination. In both processes, stabilizing selec- 
tion reduces trait divergence and diversity, and its effects on divergence are always stronger 
than on diversity. The reason will become clear in Section 4: the effective strength of selec- 
tion on trait divergence is stronger than on diversity, albeit for different reasons in sexual 
and in asexual populations. In particular, we show in Section 6 that the equilibrium ratio of 
trait divergence and diversity shows a nearly universal behavior: it decreases with increasing 
strength of selection in a predictable way, but it depends only weakly on the number of 
constitutive sites, their selection coefficients, and the recombination rate. Hence, this ratio 
provides a new, quantitative test for stabilizing selection on quantitative traits, which does 
not require genomic data and is applicable at arbitrary levels of recombination. The agenda 
of the paper is summarized in Fig. 1. For the reader not interested in any technical details, 
the summary of genome evolution (Section 2.4) together with the basics of trait statistics 
(Section 3.1) and stabihzing selection (first part of Section 4) provide a fast track to the 
selection test in Section 6. 

2 Genome evolution 

In this section, we review the sequence evolution models underlying our analysis of quan- 
titative traits. All of these models are probabilistic. They describe the dynamics of an 
ensemble of populations, any one of which is described by the frequencies of its genotypes. 
The generic genotype frequency ensemble is quite intricate, because it is neither observable 
nor computationally accessible. However, this ensemble will serve as the basis for our theory 
of quantitative traits for non-recombining traits. The ensemble description simplifies in two 
well-known limit cases: the weak-mutation regime, where a population reduces to a single 



4 




Figure 1: Evolution of a quantitative trait under stabilizing selection (schematic). The 

trait E evolves in a fitness landscape f{E) favoring a single trait value E* (red line, upper panel), 
which can be compared to neutral evolution (lower panel). Its dynamics is a stochastic process, 
which results from the underlying genome evolution (Section 2). This process can be described by 
an ensemble of populations (Section 3). An individual population from the ensemble has a trait 
distribution with mean F and variance (diversity) A; two such populations are shown as brown 
curves. The trait mean is at a distance F = F — (F) from the ensemble average (F) and at a distance 
A = T — E* from the optimal trait value E* . The phenotypic population ensemble is characterized 
by the average divergence between populations (which equals twice the ensemble variance (F^)), the 
average diversity (A), and the average distance from the trait optimum, (A). Stabilizing selection 
reduces all of these quantities compared to neutrality, but the relative change is larger for (F^) and 
(A) than for (A) (Section 4). Our theory describes a number of important characteristics of the 
population ensemble: genetic load, free fitness and predictability of evolution (Section 5). The ratio 
between divergence and diversity is the basis of a new test for stabilizing selection on quantitative 
traits (Section 6). 
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fixed genotype, and the strong-recombination regime, where genotype frequencies can be 
expressed by allele frequencies at individual genomic sites. 



2.1 Evolution of genotypes 

At the most fundamental genomic level, a population is a set of genotypes. A genotype is a 
sequence a = (oi, . . . , a^) of length I from a /c-letter alphabet (with A; = 4 in actual genomes 
and A; = 2 in our simplified models); there are K = such genotypes. In a given population, 
each genotype has a frequency > with the constraint "^^.x^ = 1. We describe the 
population state by recording the linearly independent frequencies x — {x^, . . . , x^~^) for a 
set ^ of X — 1 genotypes (the one remaining, arbitrarily chosen reference genotype slk has 
the frequency x^ = 1 — ^ae.A^'^)' 

The evolution of genotypes is a stochastic process, which generates a probability dis- 
tribution of genotype frequencies, P{x,t). This distribution describes an ensemble of inde- 
pendently evolving "replicate" populations and follows a generalized Kimura diffusion equa- 
tion [19,31], 



a.heA 



P{x,t). (1) 



Here and below, we adopt the convention that differential operators act on all functions to 
their right. The first term on the right hand side of eq. (1) accounts for stochastic changes 
of genotype frequencies by reproductive fiuctuations in a finite population (i.e., by genetic 
drift). This term is proportional to the inverse of the effective population size N and to the 
diffusion coefficients 

g-^(x) = I 'l^'^l (2) 

^ ^ 1^ a;'^(l - x*^) ff a = b. ^ ^ 

The second term describes deterministic frequency changes by mutations. In asexually repro- 
ducing populations, the coefficients m^{x) are given in terms of the mutation rates /i^ = /Xa^b 
between genotypes, 

m-(a;) = 5^(/i^xb-/x^a;-); (3) 

b 

in sexual populations, there are additional contributions from recombination. The third term 
describes natural selection. Its coefficients are fitness differences, Sb = /(b) — /{rk), where 
/(b) denotes the reproduction rate (Malthusian fitness) of a genotype b e ^ and fisLx) is 
the corresponding rate for the reference genotype slk- These k — 1 selection coefficients char- 
acterize the dynamics of the linearly independent genotype frequencies x — {x^, . . . , x^~^). 
Here we consider the simplest case, where all reproduction rates are frequency- and time- 
independent constants. In that case, the selection coefficients Sb can be written as the 
gradient of a scalar fitness landscape F{x), 

«b(x) = ^Fi^), (4) 
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which is simply the mean population fitness, 

F(a;)=/»^5]/(a)a;^ (5) 

a 

(see ref. [41] for a discussion of more general cases). Although the probability distribution 
P{x, t) of genotype frequencies gives a complete description of an evolving population en- 
semble, it is not an observable quantity. Even for moderate genome length i, there are vastly 
more possible genotype distributions x than can be recorded from the history of a single pop- 
ulation or even from an ensemble of independently evolving populations. Like the probability 
distribution over phase space in statistical mechanics, this distribution should be regarded as 
a conceptual and computational intermediate: P{x, t) is calculated using maximum-entropy 
postulates, and it is used to define and predict expectation values of observable quantities. 

Importantly, the definition of expectation values involves averaging at two distinct lev- 
els. In a given population, the genotype frequencies x determine the allele frequencies at 
individual genomic sites, 

a 

the haplotype (allele combination) frequencies at pairs of sites, 

a 

and so on, which are conveniently represented as averages of the "spin" variables 

^ f 1 if = a, 
* [ otherwise. ^ ^ 

These averages within a population are denoted by overbars. Connected correlation functions 
at a single site, 

K = (<^f-y?)«-?/f) = (1 - yt), (9) 

are components of the sequence diversity TTj — Ylt.=i ^f^ correlations between different sites, 

^S^«-yfM-y') = yS-yty^ i^^j)^ (lo) 

measure linkage disequilibrium, i.e., biases in the association of alleles to haplotypes within 
a population. For all of these quantities, the genotype frequency distribution P{x, t) defines 
expectation values in an ensemble of independently evolving populations, 

( <^ > = j afa'j... P{x, t) dx; (11) 

averages across populations are denoted by angular brackets, (■). Such nested correlation 
functions can often be decomposed into independent fiuctuation components within and 
across populations; for example, 

< ^ > - ( ^ ^ > + < « - yt) (^^ - y'j) > - {yty') + ■ (12) 
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In particular, fitness interactions (epistasis) can generate allele frequency correlations {Viy]) 
even if linkage disequilibrium vanishes. 

All frequency correlation functions of the form (11) can, in principle, be computed from 
the solution of the diffusion equation (1). This is impossible in practice, however, because no 
general analytical solution exists. In particular, the distribution P{x, t) does not converge to 
an evolutionary equihbrium, which is defined as a state with detailed balance (see ref. [40] 
for a review of detailed balance in an evolutionary context). We remind the reader that a 
diffusion equation of the form (1) has an equilibrium distribution if and only if the vector 
field v^{x) = Sbe^ {rrv^{x) + g^^{x)s\i{x)) satisfies the integrability conditions, 

d , d 

Q^^9e.a'{x)v^'{x) - Q^^9he.'{x)v^\x) = 0, (13) 
a' a' 

which imphes that v°'{x) can be written as the gradient of a scalar function. It is easy to see 
that the frequency-dependence of the diffusion matrix g^^{x) makes already the mutation 
vector field m^(a;) non-integrable. Hence, even if the selection coefficients s-^{x) are the 
gradient of a scalar fitness landscape as given by eq. (4), there is no general evolutionary 
equilibrium. We now discuss the two known special cases in which the diffusion equation 
(1) does have a solvable equilibrium, which is the analogue of the Boltzmann equilibrium in 
statistical thermodynamics. 

2.2 Weak- mutation regime 

This regime is defined by a low genome- and population-wide mutation rate per generation, 
liNl ^ 1 [25]. With typical values 1.1N ~ 10^^, this regime applies to short sequence 
segments with a length up to about 10 base pairs, which are the genomic basis of microscopic 
traits. Transcription factor binding sites in prokaryotes and simple eukaryotes with a typical 
length of about 10 base pairs are examples of this kind [38,39]. In such segments, a single 
genotype is fixed in the population at most times. This genotype evolves through occasional 
polymorphisms at a single genomic site, but co-occurrence of polymorphisms at multiple 
sites can be neglected. We can then project the Kimura equation (1) onto a Master equation 
on the space of fixed genotypes, 

^P(a,i) = J][Mb^a^(b,i) -iia->b^(a,t)]. (14) 

b 

The substitution rates are given by the classic Kimura-Ohta formula [30,32], 

«a^b = /Xa^b2iV(/(b) - /(a))/(l - exp[2iV(/(b) ~ /(a))], (15) 

with selection coefficients given by the discrete fitness landscape /(a). We make an assump- 
tion on neutral evolution: it occurs by point mutations with site-independent rates iia^hi 
which satisfy the detailed-balance relations 

Po{a) iia^b = Po{b) Hb-^a- (16) 
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This detailed-balance assumption, which is part of all standard neutral mutation models, 

reduces the number of independent rate constants from 12 to 9. The resulting equilibrium 
single-nucleotide distribution p{a) (a = 1, . . . ,k) describes the effect of mutational biases 
(if all rates are symmetric, Ha^b = A*6^a, it leads to a fiat single-nucleotide equilibrium 
Po{o,) = ^/k). The detailed-balance condition in the weak- mutation regime is much weaker 
than the corresponding condition (13) for frequency evolution, which constrains an entire 
function v'^{x). 

In an arbitrary fitness landscape /(a), the full dynamics (14) with (16) has an equilibrium 
probability distribution of fixed genotypes [7, 49] 

Peq(a) = ^Po(a)exp[27V/(a)], (17) 

which is the product of the factorizable neutral equilibrium 

e 

Po(a) = l[po(ai) (18) 

i=l 

and the Boltzmann factor exp[2A'^/(a)]. Here and below, Z denotes a normalization factor. A 
generic fitness landscape generates cross-population allele correlations {yiVj) between sites, 
but linkage disequilibrium vanishes without any assumptions on the recombination rate. 



2.3 Strong-recombination regime 

In this regime, linkage correlations become small because of rapid allelic reassortments in 
the population. We can then approximate the frequency of a genotype by the product of its 
allele frequencies, 

x^^yr-.-yT- (19) 
This approximation, which we will refer to as free recombination, describes complete link- 
age equilibrium. It becomes exact in the limit of infinite recombination rate and infinite 
population size. Linkage equilibrium is the standard assumption of quantitative genet- 
ics [20,22,33-35]; it is often applied to large genomes in sexually reproducing populations, 
which are the genomic basis of macroscopic quantitative traits. Given the factorization (19), 
we can project the Kimura equation (1) onto a diffusion equation for the joint distribution of 
allele frequencies. In the simplest case of a two-letter genomic alphabet, this equation takes 
the form 

° 1 d"^ d 

i9{yi) - -K-{m{yi) + g{yi)si{y)) 



dt 



2Ndyr'''' dy, 



P{y,t). (20) 



Here, y = {yi, . . . ,ye) denotes the set of allele frequencies, g{y) = y{l — y) and m{y) = 
li{l — 2y) are the diffusion and mutation coefficients, and Siiy) = dF{y)/dyi are the selection 
coefficients for alleles, with F{y) = ■■■yT- arbitrary fitness landscape 

F(y), the projected Kimura equation has an equilibrium distribution [54,55], 

Peq(y) = ^Po(z/)exp[2iVF(|/)], (21) 
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which is the product of the factorizable neutral equihbrium 

Po{y)-^U[y^{l-y^r (22) 

and the Boltzmann factor, exp[2A'^F(y)]. In this case, equihbrium emerges because the 
neutral distribution is the product of one- dimensional allele frequency distributions, for which 
the integrability condition (13) is always fulfilled. Just as in the weak-mutation regime, a 
generic fitness landscape generates allele frequency correlations {yfi/j), which are compatible 
with linkage equilibrium. The strong-recombination calculus can be extended to populations 
with a large but finite recombination rate [17,42]. Such populations still reach an evolutionary 
equilibrium of the form (21); however, the neutral distribution Po{y) no longer factorizes and 
there are small but systematic linkage correlations (Tr^j) [42]. 

2.4 Summary 

Genomic evolution under mutations, recombination, genetic drift, and selection can be de- 
scribed by a Kimura diffusion equation at the genotype level [19,31]. The general Kimura 
equation does not have a closed solution. In the regimes of low mutation rate or high recom- 
bination rate, where linkage disequilibrium is small, we can project this equation onto fixed 
genotypes or allele frequencies, respectively. These projections are shown in Table 1. The 
projected equations have solvable equilibria of the form P = Pq exp[2N F]; sec equations (17) 
and (21). The "Boltzmann" factor exp[2A^F] links the equilibrium probability distribution 
under time-independent selection, P, with the corresponding distribution for neutral evolu- 
tion, Pq] this relation can serve as a starting point for the inference of selection. However, 
genomic equilibria do not exist for strongly coupled multi-site evolution with large linkage 
disequilibrium, which is common in asexual populations and even in sexual populations for 
compact, intermediate-size genomic regions [12]. 

3 Evolution of quantitative traits 

In this Section, we first introduce the basic statistical observables for quantitative traits. In 
the low-mutation and in the strong-recombination regime, we obtain phenotypic equilibria by 
projection from the genomic equilibrium distributions discussed in the previous Section. For 
complex traits evolving under linkage disequilibrium, we show that projection of the general 
genotype dynamics leads to a stationary non-equilibrium distribution, which describes the 
statistics of trait divergence and diversity under time-independent selection. Further projec- 
tions lead to equilibrium marginal distributions of trait divergence and diversity, which are 
the basis for our subsequent analysis of stabihzing selection. The projections from genomic 
to phenotypic distributions are also shown in Table 1. 
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allele frequencies 
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diversity 
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q{£) 


P.4XA) 
/ \ 

Pc4T) Pcq(A) 


P.UTA) 
/ \ 

Pc4T)P,^(^) 



Table 1: From genotypes to phenotypes. This table shows the genomic and phenotypic 
stationary population ensembles discussed in the text. These ensembles are obtained by differ- 
ent projections of the diffusive genotype dynamics (1) under time-independent selection, which are 
marked by arrows. In the low-mutation regime, we obtain an equilibrium distribution of fixed 
genotypes, which can be projected further onto an equilibrium of fixed trait values; see equations 
(17) and (33). In the strong-recombination regime, we obtain an equilibrium distribution of allele 
frequencies, which can be projected further onto a joint equilibrium of trait mean and diversity; see 
equations (21) and (34). For complex traits evolving without recombination, we obtain a station- 
ary non-equilibrium distribution of trait mean and diversity, which can be projected further onto 
equilibrium marginal distributions; see equations (42), (48) and (52). 

3.1 Trait statistics within and across populations 

The subject of this paper is the evolution of quantitative traits with a heritable component, 
which depends on an individual's genotype. Here we study the simplest case of an additive 
map from genotype to phenotype, and we assume a binary genomic alphabet (extension to 
a /c-letter alphabet is straightforward). Any phenotype E can then be written in the form 

E(.)^±E,a, w,th...{j l;;^ (23) 

i=l ^ 

Here the phenotype is measured from its minimum value, a = (ai,. . . ,0^) is the genomic 
sequence at its constituent sites, a* is the allele conferring the larger phenotype at a given 
site, and E'j > is the phenotypic effect at that site, i.e., the difference in trait value between 
the two alleles. We define the allelic average Fq and the overall effect amplitude Eq by 

To-^E^- ^o-^E^^f- (24) 

i=l i=l 

We are interested in the evolution of complex molecular traits, which depend on multiple 
genomic sites. If the number of constituent sites is sufficiently high (such that is of order 
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1 or larger), such traits are generically polymorphic in a population, even if most individual 

sites are monomorphic (i.e., 9 = jiN ^ 1, which is the case in most populations). The trait 
values in the individuals of a population follow a distribution W(-E'). Here we parametrize 
this distribution by its mean and its variance, which is called the trait diversity [4, 10, 17,34]: 

V = E^^E{ai)x^, A = E^-V'' ^^Elx^-^E^Ey,x^x^ (25) 

a a a,b 

Using eq. (6), the trait mean can be written as a function of the allele frequencies, 

i £ 

V{y) = Y,Eia, = Y,Eiyi (26) 

i=l i=l 

The trait diversity can be decomposed into the additive trait diversity Ai{y), which depends 
only on the allele frequencies, and the trait autocorrelation A2(7r), which depends on the 
linkage disequilibria between the constituent loci, 

e 

I 

= J]e2|/,(1-|/,) + J]E,E^.7r,, = Ai(|/) + A2(7r), (27) 

where we have used eqs. (9) and (10). As will become clear, the trait diversity is therefore 
more strongly affected by linkage and recombination than the trait mean. In the strong- 
recombination approximation of quantitative genetics, the population is at linkage equilib- 
rium and the trait diversity reduces to its additive part, A ~ Ai. Under finite recombina- 
tion, stabilizing selection generates a negative trait autocorrelation A2, and the assumption 
of linkage equilibrium will lead to an overestimation of A. 

Similarly to genotype evolution, the stochastic evolution of a quantitative trait generates 
a probability distribution (5(r, A, t), which describes an ensemble of independently evolving 
populations, each having a trait distribution with mean F and variance A; see also [10]. The 
probability (5(r, A, t) is a sum of probabilities of genotype frequencies, 

g(r,A,i) = y h{Y{y{x))-Y)h{A{y{x),i:{x))-A)P{x,t)dx, (28) 

where 5(-) is the Dirac delta function. Furthermore, we assume that selection acts on a 
trait's constituent genotypes only via the trait itself, 

/(a) = /(£;(a)); (29) 

that is, all genotypes with the same trait value E have the same fitness /(-B). The genotypic 
fitness landscape F{x) given by eq. (5) then defines a phenotypic fitness landscape 

F(F, A) ^ /(F, A) = /(F) + ^Ar(F), (30) 



12 



which contains the leading terms in the Taylor expansion of f{E) around the trait mean 

r. The phenotypic evolutionary scenario is illustrated in Fig. 1 for evolution in a single- 
peak landscape f{E) and for neutral evolution (these cases are analyzed in detail in the 
next section). Each population drawn from the ensemble distribution Q{r,A) has a trait 
distribution with mean F and variance A. The trait mean of a given population is at a 
distance 

f = F - (F) (31) 
from the ensemble average (F) , at a distance 

A = T-E* (32) 

from the optimal trait value E*, and two populations have a square trait distance (Fi — F2)^, 
which is called their trait divergence. Statistical theory predicts ensemble averages such as 
(F^), (A), (A), and ((Fi — F2)^) = 2(r^). In this paper, we focus on stationary ensembles 
under time-independent selection. The trait divergence can also be defined for a single 
population at two different times, D(t2 — ti) = (F(t2) — F(ti))^, or more generally for two 
populations with a common ancestor. Under time-independent selection, D[t) reaches the 
equilibrium ensemble divergence for long times, limf_!.oo(-D(t)) = ((Fi — F2)^). The statistics 
of time-dependent trait divergence will be analyzed in another paper [28]. 

The projection from genotypes to phenotypes given by eqs. (28) and (30) can immediately 
be put to use in the regimes of low linkage disequilibrium discussed in the previous section, 
where an evolutionary equilibrium exists at the genomic level. In the weak-mutation regime, 
we obtain an equilibrium distribution of fixed phenotype values by projection from eq. (17), 

Q,^{E) = 1 QoiE) exp[2Ar/(E)], (33) 

this type of equilibrium distribution has been used in refs. [7,38,39]. In the strong-recombination 
regime, the phenotypic equilibrium obtained by projection from eq. (21), 

Qeq(r, A) = ^go(r, A)exp[27VF(F,A)]. (34) 

has been analyzed in detail in refs. [4, 17]. 

3.2 Joint evolution of trait mean and diversity 

As discussed in the previous section, this equilibrium calculus is not applicable to correlated 
evolutionary processes in non-reconibining or slowly recombining genomes, which evolve 
large values of linkage disequilibrium. To analyze such processes, we proceed differently: 
we directly use the Kimura equation for genotypes to obtain by projection a self-consistent, 
approximate diffusion equation for the phenotypic ensemble distribution Q{r,A,t). In this 
paper, we study the case of strictly asexual, non-recombining populations. By projection 
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from eq. (1), we find the phenotypic diffusion equation 



d_ 
di 



Q{r,A,t) = 



2N \ dr^ 

d_ 

~dK 



9 



d_ 



/ A I AA \ 

[m +g sa) 



Q{r,A,t) 



(35) 



with diffusion coefficients g^^\g^^, mutation coefficients m^, m^, and selection coefficients 
sr,SA that depend on the variables T and A. We obtain the diagonal diffusion coefficients 



9 



rr 



^ dr dr 



a,b 



9 



E{a)E{h) [-x^x^{l - 5l) + x^{l - x^)5l] 

a,b 

{E - r)2 = A, 



(36) 



,AA 



^ OA OA ,b 



a,b 



dx^ dx^ 



9 



= Y^i^i^T - '2EE{si)){El - 2EE{h)) [-x"x^(l - S^) + x^(l - x^)5^] 

a,b 

= {E - r)4 - A^ fti 2A2. 



(37) 



These diffusion coefficients reflect stochastic changes in trait mean and diversity by sampling. 
It is clear that the fluctuation amplitude (36) for the trait mean is set by the trait diversity. 
The corresponding amplitude (37) for the trait diversity is specific to asexual evolution: 
sampling of a set of complete genotypes with trait values E^ from a Gaussian distribution 
yV{E) with variance A leads to a distribution of sample variances with variance 2A'^. This 
relation changes in recombining populations, where sampling is broken down to individual 
alleles. For more general trait distributions W{E), the amplitude g'^'^ given by eq. (37) 
involves higher moments [17,42]; that is, the closed form (1) of the dynamics for F and 
A is a truncation. As shown by our numerical results, this truncation leads to accurate 
approximations for complex quantitative traits, because their actual trait distribution is 
approximately Gaussian. As we have anticipated in writing eq. (35), off-diagonal diffusion 
can be neglected by symmetry. 



.rA 



dA dr 
dx^ dx^ 



ab 



= {E^(a) - 2EE{si))E{h) [-x^x^{l - S^) + x^{l - x^)5^] 
= {E - F)3 ^ 0. 



(38) 
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This coefficient would lead to additional terms, such as {2N)-^{d'^ /drdA)g^^Q{r, A,t). 
The mutation coefficients are 



r 

m 



= J] -2y,) = -2/x(r- To) (39) 

i=l 

= -MA-E',)-^ + 0{e'), (40) 

where Fq and Eq are given by eq. (24). The term A/A?" in (40) appears due to the nonlinear 
dependence of the trait diversity on the allele frequencies yi (see, e.g.. Chapter 4 of ref. [23]). 
Finally, the selection coefficients are the gradient of the phenotypic fitness landscape (30), 

«r = ^F(r,A), .A = ^F(F,A). (41) 

The two-dimensional diffusion equation (35) gives a closed, analytical description of trait 
evolution under complete genetic linkage. As we show in the next section, it provides nu- 
merically accurate results at least for the marginal distributions Q{T) and Q(A) over a wide 
range of evolutionary parameters. However, it has the same basic difficulty as the geno- 
typic Kimura equation (1): it does not have an equilibrium solution, because the mutation 
coefficient field is non-integrable, d{{g^^)~^mF) / dA — d{{g^^)~^m^) / dV ^ 0. In the Ap- 
pendix, we show that equation (35) leads instead to a non-equilibrium stationary distribution 
Qstat(r, A), which is shown in Fig. 2. This distribution satisfies the scaling relation 

g,tat(r,A) = r3/2 4,,,(r^/2(r- (r)),r^A) (42) 

for large values of ^, with ensemble averages (F) and (A) of order I. According to this 
relation, the average (F) and the fluctuations F = F — (F) of the trait mean in the stationary 
ensemble scale in accordance with the central limit theorem, 

(F)~£, (f")~r/' (n = 2,3,...), (43) 

which implies that fluctuations become subleading in the large-£ limit, F = (F) ± 0(£^/^). 
This scaling also occurs in sexual populations. It is analogous to the thermodynamic limit 
for macroscopic systems, which is familiar in statistical thermodynamics [17]. However, the 
average (A) and the fluctuations A = A — (A) of the trait diversity scale in a different way, 

(A) - £, (A") ~ r (n = 2, 3, . . . ). (44) 
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Figure 2: Non-equilibrium stationary trait distribution under complete genetic link- 
age. Stationary joint distribution of trait mean and diversity, Qsta.t(r, for a non-recombining 
population in a quadratic fitness landscape. The figure shows simulation results for a quantita- 
tive trait with £ = 100 constituent sites of equal effect. The distribution Qsta,t(^,^) is Gaussian 
in the T direction, but strongly non-Gaussian in the A direction (the resulting marginal distribu- 
tions are shown in Fig. 3). It maintains a stationary probability current, which is shown in Fig. 7. 
Other system parameters: neutral sequence diversity 9 = = 0.0125, scaled fitness landscape 
2Nf{E) = c{E - E*f/El of strength c = 2.5 with a fitness optimum E* = Q.^Tq. 

This scale-invariance of the trait diversity statistics in large-^ limit is a consequence of co- 
herent, genome-wide linkage disequilibrium fluctuations in the absence of recombination. It 
is generated by sampling from a set of genotypes with trait values Eg, from a distribution 
W(-E) with variance A ~ There is no central-limit theorem, because the number of 
these genotypes grows only weakly with ^ [50]. In contrast, fast recombination generates a 
number of genotypes that grows exponentially with £, which leads to the standard scaling 
(A") ~ given by the central limit theorem (see [17] and the discussion in the next sec- 
tion). These differences in fluctuation statistics are mirrored by the properties of population 
genealogies: for asexual evolution, there is a single genome-wide genealogy of all genotypes. 
Standard coalescence theory then predicts diversity fluctuations distributed exponentially, 
with variance proportional to the square of the coalescence time, which of order A^^, and 
to the square of the genome- wide mutation rate, which in turn is proportional to In 
contrast, recombination generates many parallel genealogies, which average out the diversity 
fluctuations. 

Because the joint evolution of trait mean and diversity is a non-equilibrium process, the 
diffusion equation (35) does not have a simple analytical solution. We now project this 
dynamics further onto its marginals for T and A. The ensemble distributions Q{T,t) and 
Q(A, t) follow coupled one- dimensional diffusion equations, which turn out to have analytical 
equilibrium solutions. 

3.3 Marginal evolution of the trait mean 

By integrating over the trait diversity in eq. (35), we obtain a one-dimensional diffusion 
equation for the trait mean. This integration amounts to replacing the variable A, which 
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appears in the diffusion coefficient g^^ and in tfie selection coefficient sp, by its expectation 
value (A) . The projected equation reads 



d_ 
di 



Q{r,t) 



9 



rr 



d 



2N dr^ dr 



Q{T,t) 



with the effective diffusion coefficient 



(45) 



(46) 



the mutation coefficient = — 2/i(r — (r)o) given by eq. (39), and the selection coefficient 
sr, which is the gradient of the effective fitness landscape 



F(r) = /(r,(A)) = /(r) + -(A)r(r). 



This equation has an equilibrium solution 



1 ^ 



Qeq(r) = -go(r)exp[2iVF(r)], 



with 



Qo(r) 



29 
ttTA) 



exp 



1 (r - Tc 



2 (A)/4^ 



(47) 



(48) 



(49) 



and To given by eq. (24). The Gaussian form of Qo(^) is valid for sufficiently large values 
of £. It implies the scaling form (43) of average and fluctuations of F, in accordance with 
the central limit theorem. Since F depends only on the allele frequencies of the constituent 
loci and not on their linkage correlations, the diffusion equation (45) and the form of its 
solution (48), (49) are valid regardless of recombination. However, the distribution Qo{T) 
depends on the average diversity (A) under selection, which enters the effective diffusion 
coefficient (46). Hence, Qo(r) differs from the neutral distribution QqIT). Because (A) 
depends on recombination (see below), the statistics of the trait mean also acquires a small 
but systematic dependence on the recombination rate. 



3.4 Marginal evolution of the trait diversity 

For non-recombining populations, we obtain a one- dimensional diffusion equation for the 
trait diversity from eq. (35), 



d_ 
di 



Q{A,t) 



1 

27V aA2 



,AA 



9 



d_ 
dA 



( A I AA~ \ 

[m +g sa) 



Q{A,t) 



(50) 



m 



with the diffusion coefficient g^^ = 2A^ given by (37), the mutation coefficient 
— 4/^(A — Eq) — A/N given by (40), and the selection coefficient sa, which is the gradient of 
the effective fitness landscape 

F{A) = l{f"{T))A. (51) 
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This equation has an equihbrium solution 



Qeq(A) ^ ^go(A)exp[27VF(A)], (52) 



where Qq{A) is the neutral equilibrium 



^0 



4:9E, 



(no recombination). (53) 



with the normahzation Zo — {29 E^) ^ rEuier(2 + W). This distribution has mean and 
variance 

(A)o = AeEl{i-Ae) + o{e^), 

((A-(A)o)^)o = AeEl{l-W) + 0{e^) (no recombination). (54) 

It is of the form (5o(A) = ^~^Qo{^~^A) with a scale-invariant shape function Qo, which 
imphes the coherent scahng (44) of diversity mean and fluctuations (see also Appendix) . 

The trait diversity equilibrium (52, 53) can be compared with its counterpart for free 
recombination. The equilibrium distribution Qcq(A) for the free-recombining traits is also 
of the form (52), with the neutral distribution (5o(A) obtained by projection from the allele 
frequency distribution (22), 

go(A) = b{A-J2EM^-yi))Po{y)dy^...dy, 

i=l 

(free recombination). (55) 
For sufficiently large i, this distribution is again Gaussian with mean and variance [17], 
(A)o = A0E^^{1-Ae) + O{0'), 

((A-(A)o)^)o = ^ - ^) +C'(^^) (free recombination), (56) 

i=l ^ ' 

which implies the standard scaling of diversity average and fluctuations, (A) ~ i and (A") ~ 
^"/2for (n = 2,3,...). 

In a generic fitness landscape, the equilibrium distributions (48) and (52) for trait mean 
and diversity depend on each other, and consistent joint solution has to be obtained itera- 
tively. Mean and diversity decouple in a linear fitness landscape [4], and the dynamics of the 
diversity is still autonomous in a quadratic fitness landscape. This case will be discussed in 
the next section. 

We test our analytical results by simulations of a Fisher- Wright process under stabi- 
lizing selection and at neutrality. We evolve a population of N individuals with genomes 
a.^^\ . . . , a^^\ which are bi-allelic sequences of length i. A genotype a defines a phenotype 
E{a.) — Yli=i Eiai] the phenotypic effects Ei drawn from various distributions. In each 
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generation, the sequences undergo point mutations with a rate t/i per generation (where r 

is the generation time). The sequences of next generation are then obtained by multino- 
mial sampling; the sampling probability is proportional to [1 + r/(i?(a)] with the fitness 
f{E) = —Co {E — E*Y (details are given in the next section). For sexual populations, we 
permute the alleles Oj^i, . . . , Cj^jv at each genomic site i between the individuals in each gen- 
eration, which amounts to recombination with an infinite rate. As shown in Fig. 3, the 
analytical equilibrium distributions Qeq(r) and Qq{T) given by eqs. (48, 49), as well as the 
distributions (5cq(A) and Qo(A) given by eqs. (53, 55), are in good agreement with simulation 
results. Stabilizing selection shifts the average and reduces the variance of the distribution 
(5(r), and it reduces average and variance of the distribution (5(A) compared to neutral 
evolution. The dependence of these effects on the strength of selection is analyzed in the 
next section. 

4 Trait equilibria under stabilizing selection 

We now apply our statistical model to quantitative traits under stabilizing selection, a sce- 
nario described by evolutionary equilibrium in a quadratic fitness landscape, 

f{E) = r-co{E-E*f, (co>0). (57) 

This scenario is probably a reasonable approximation for many actual traits, which have 
high-fitness values in a certain range around their optimum value E* [1,4,17]. For example, 
it applies to the expression level of a gene: small changes in expression may be buffered 
by compensatory changes in the regulatory network and will affect fitness only weakly, but 
larger changes are often deleterious, as it is evident from the large number of genetic disorders 
associated with gene copy number variation. 

Stabilizing selection changes the distribution of trait values in a population, VV(-E), which 
can be parametrized by changes in the trait mean F and the diversity A. Statistical theory 
describes the expectation values of these changes in an ensemble of populations. At a quali- 
tative level, the main effects are already clear from the previous section: stabilizing selection 
decreases the average squared distance from the fitness optimum, (A)^ = ((F) — E*Y, the 
average equilibrium divergence between populations, ((Fi — F2)^) = 2(F^), and the average 
diversity, (A). We now derive analytical expressions for these effects in non-recombining 
populations and under free recombination, and we analyze their dependence on the selection 
strength cq (the fitness maximum /* is an arbitrary constant, because the evolution equation 
(35) depends only on fitness gradients); see also rcfs. [1,9,17,56,57] for effect of stabilizing 
selection on free-recombining macroscopic traits. Compared to a generic fitness landscape, 
the analysis is somewhat simplified for a quadratic fitness landscape (57), because the mean 
population fitness separates, 

F(F, A) = r - co(F - E*f - coA. (58) 

For a quantitative analysis, it is useful to measure phenotypes in a natural unit, which 
avoids the arbitrariness of fixed units (such as centimeters or inches for body height). Here 
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Figure 3: Equilibrium trait distributions under stabilizing selection and at neutrality. 

(a) Equilibrium distribution (5eq(r) of the trait mean in a quadratic fitness landscape (filled circles) 
and corresponding neutral equilibrium Qo(r) (empty circles) (green: no recombination, blue: free 
recombination), (b) Equilibrium distribution (5cq(A)of the trait diversity (green: no recombination, 
blue: free recombination). Theory predictions for these distributions are shown as solid and dashed 
lines. System parameters: £ = 100 trait loci of equal effect Ei = \ {i = 1 . . . t)^ neutral sequence 
diversity 9 = = 0.0125, scaled fitness landscape 2Nf{E) = -c{E - E* f/El of strength c = 5 
with a fitness optimum E* = 0.7L. Results for other effect distributions are shown in Fig. 6. 
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we express trait values in units based on the effect amplitude (24), 

E r , A 



e=^, ^=^' '^=^' ^^^) 

-C/n -C/n -C/n 



and in the same way e* = E*/Eq, A = A/£'o and 7 = 'T/Eq. These scaled values are pure 
numbers (we distinguish them by use of lower case letters from the raw data). The scaling 
(59) has a straightforward biological interpretation: Eq is the trait variance in an ensemble of 
random genotypes, which would result from neutral evolution in the weak-mutation regime, 

£;o^ = iim((r-(r))V (60) 

We also use the effect amplitude to define the scaled strength of stabilizing selection, 

c = 2NEl Co, (61) 

which can be interpreted as the difference between the fitness maximum /* and the average 
fitness in the random ensemble, 

c = 27V/* - hm (27V/)o , (62) 

where fitness (growth rate) is measured per 2N generations and we have assumed that 
selection does not shift the trait average (i.e., E* = (r)o). Such fitness differences are 
referred to as genetic load, which is discussed in section 5.1. 

4.1 Trait average under stabilizing selection 

In a quadratic fitness landscape, the equilibrium distribution of the trait mean is Gaussian 
for sufficiently large values of i, 



Qeq(r) = ^Qo(r)exp[27V/(r)] = ^exp 



Zr Zr l{d) 



(63) 



with 7o = X]i=i^i/2! as given by eqs. (47), (48), and (49) and Z^ as the appropriate nor- 
malization factor. This distribution has the scaled moments 

(f>.(,V<T>^ ^ (f).^^^^^. (64) 

with (A)o = 7o — e* and (7^)0 = 1 — 4^^ + 0{d'^). In the regime of weak selection (c <C 1), 
these moments depend on the selection strength c in a universal way, 

(A)2 = (A)^(l-4c) + 0(c^c/£,c^), 

(f) = (f)o(l-2c) + 0(c^c/£,c^), (65) 
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Figure 4: Trait moments under stabilizing selection, (a) The squared average distance of 
the trait mean from the fitness optimum, (A)^, (b) the variance of the trait mean, (7^), which equals 
half the average equilibrium divergence, and (c) the average diversity (6) are plotted against the 
selection strength c (green: no recombination, blue: free recombination). Other system parameters 
are as in Fig. 3. All quantities are scaled by the effect amplitude Eq. In both recombination regimes, 
the effect of stabilizing selection on the trait diversity is seen to be smaller than on the trait mean. 



22 



because the effect of selection on the trait diversity is subleading {(5) / {S)o = 1 + 0{c/£, c9), 
see eq. (72) below). For larger values of c, these moments acquire a noticeable dependence 
on (6), and thereby on the recombination rate. For asexual populations, we obtain the 
strong-selection regime {c9 ':$> 1) 

if) = 7^ [1 + c-^/2)] (no recombination), (66) 

2c 

where we have used eq. (68) below. Evaluating this regime does not make sense in the 
free-recombination approximation, because if epistatic selection is strong, the assumption of 
linkage equilibrium breaks down for any finite recombination rate. 



4.2 Trait diversity under stabilizing selection 

The equilibrium distribution of the trait diversity is 



geq(A) = — go(A)exp(-27VcoA) 



(67) 



with (5o(A) given by (53) and (55) and Z/^ as the appropriate normalization constant. This 
distribution does not depend on the statistics of the trait mean and determines the scaled 
average diversity in asexual populations by. 



{5) 



1 



r2-4^xp(-^-c5) d5 




W A;i+49[4V^c] 



= {6)o[i + g{ec)] 



(no recombination). 



(68) 



where kn{z) denotes the modified Bessel function of the second kind. The average diversity 
of free-recombining traits reads (see also ref. [17]), 



2\2e 



exp(— c(5j) d5i 



Fiji + 20,3/2 + 26,-0 ef/4] 
Fi[2^,l/2 + 2^,-c ej/A] 



9 
2 



de K,{e)e 



2 Fi [1 + 29, 3/2 + 29, - (c/i) eV4] 
F,[29,l/2 + 29,-{c/£) 6^4] 



free 



(free recombination). 



(69) 
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where Fi[a, b, z] is the regularized confluent hypergeometric function and we have introduced 
the effect density 

«:(6)^^X]6(6-£V2e^). (70) 
* i=i 

We can again expand these expressions to leading order in c, 

(S) = {S)o[l-^ec + 0{e^c^)] (no recombination), (71) 



1 _ '^^^nf- — 



(free recombination), (72) 



where /t„ denotes the n-th moment of the distribution K{e) {n = 1,2,...). For asexual 
populations, we obtain the strong-selection regime {c9 3> 1) 



(*> = (■5)0 



(no recombination). (73) 



Again, evaluating this regime does not make sense in the free-recombination approximation, 
because approximate linkage equilibrium cannot be maintained at any finite recombination 
rate. For c/i » 1, selection changes even qualitatively: it becomes balancing at individual 
trait loci and would act to increase the trait diversity. 

Our analytical results (64), (68), and (69) for trait equilibria under stabilizing selection 
are shown in Fig. 4 together with numerical simulations. As expected, the behavior of 
the trait diversity depends more strongly on the recombination rate than that of the trait 
mean. However, there is an important and universal feature: stabilizing selection affects the 
trait diversity always less than its mean. This feature, which will be the basis for a test of 
stabilizing selection on quantitative traits, is explicitly demonstrated by our solution. As 
shown by eqs. (68) and (69), selection on trait diversity has an effective strength 

9c <^c (no recombination), (74) 

c/i <^ c (free recombination), (75) 

which involves a small pref actor compared to the selection strength c acting on divergence. 
These prefactors reflect different mechanisms of stabilizing selection acting on trait diversity. 
In asexual populations, selection acts on a distribution of genotypes, which generates a 
neutral trait diversity by a factor 9 smaller than the neutral trait divergence. In sexual 
populations, selection acts on individual trait loci, and the mean square trait amplitude of 
an individual locus by a factor of order (l/i) smaller than the mean square amplitude Eq of 
the entire trait. 



5 Fitness and entropy under stabilizing selection 

The distributions of trait mean and diversity derived in the previous section also determine 
the fitness and entropy statistics in the equilibrium population ensemble. This statistics 
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provides a few biologically relevant numbers: it quantifies how well adapted typical popula- 
tions are under stabilizing selection, how much adaptation has occurred between neutrality 
and the adapted state, and how much measurements in one population can predict about 
another population evolving in the same fitness landscape. 



5.1 Genetic load 

How far away is a population from the fitness peak? This question is answered by the genetic 
load 

L = r- 1 (76) 

which is defined as the difference between the fitness maximum and the mean population 
fitness (and is conveniently measured in units of 1/2N) [13,14,26,37]. In the quadratic 
fitness landscape (57), we can decompose L into a component associated with the trait 
mean, 2NLy = c(7 — e*)^, which is generated mainly by substitutions away from the fitness 
optimum, and the diversity load, 2NLa = cd, which is generated by trait polymorphisms. 
Our statistical theory predicts the ensemble average of the genetic load at equilibrium, 

{2NL)^c{{Xf + {f) + {S)), (77) 

in terms of the leading moments of trait mean and diversity, which are given by eqs. (64), 
(68), and (69). Fig. 5(a) shows that the total load and its two components depend on the 
strength of selection in a non-monotonic way. For weak selection, the main load component is 
{2NLr), but {2NL/s) dominates for strong selection. This reflects our result that stabilizing 
selection affects the trait diversity less than its mean. 



5.2 Free fitness and fitness fiux 

How far away is a population ensemble from neutral evolution? This can be measured in 
two ways: by the difference in average scaled fitness between that ensemble and the neutral 
ensemble 

{2NJ)q - {2Nf)o = {2NL)o - {2NL)q, (78) 

and by the relative entropy or KuUback-Leibler distance between the ensemble under selec- 
tion and the neutral ensemble. 



H{Q\Qo) = J dTdAQ{T, A)log 



Q(r,A) 



(79) 



.Qo(r,A). 

The difference between scaled fitness and relative entropy is called free fitness, 

F{Q) = {2Nf)Q-H{Q\Q,)- (80) 

see refs. [3,7,29,41,49]. This quantity is of particular importance, because it satisfies a 
growth principle similar to Boltzmann's iJ-theorem in statistical physics: for any evolution- 
ary process in a time-independent fitness landscape which has an equilibrium, the free fitness 
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Figure 5: Genetic load, fitness flux, and predictability of evolution, (a) The average 
genetic load (L) (full lines) with its components (Lr) (dotted lines) and (La) (dashed lines), (b) the 
equilibrium fitness flux <I>eq with its components $eq,r (dotted lines) and $eq,A (dashed lines), and 
(c) the predictability V are plotted against the selection strength c (green: no recombination, blue: 
free recombination; fitness is measured in units of 1/2 A^). See definitions in eqs. (77), (82), and 
(85). Other system parameters are as in Fig. 3. 



J^{Q{t)) increases monotonically with time and has its maximum at equilibrium [29,41,49]. 
Here we approximate the stationary trait distribution under stabilizing selection by the prod- 
uct of its equilibrium marginal distributions, (5stat(r, A) ?a Qeq{T)Qcq{A) = Qcq, the same 
approximation is used for the neutral distribution Qo{r,A) (the results in the Appendix 
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show that this is numerically justified). We then obtain the relative entropy 

H{Qec,\Qo) = -c((A)2 + (f )) - logZr - c{6) - logZ^ (81) 
and the difi^erence in free fitness or fitness flux 

2iV$eq = -^(Qeq)-^(Qo) 

= -c((A)2 + if) + {5)) - i/(QeJQo) + c{{\)l + (f )o + (5)o) 

= c{{X)l + {f)o) +logZr + c{6)o + logZA, (82) 

with logZr ~ {X)l{e - (^c)V2) _ (1/2) lege and log^A ~ -4(0c)i/2 + (3/4) log c. The 
scaled fitness flux 2A^$cq measures the total amount of adaptation between the neutral 
equilibrium and the equilibrium under stabilizing selection^ [41]. As shown in Fig. 5(b), this 
fiux is always positive and increases with the selection strength c. Similarly to the genetic 
load, it can be decomposed into contributions of the trait mean and the trait diversity, 
2N^ = 2iV$eq,r + 2A/"$eq,A. The term 2A^<l>cq,r = c((A)g + (7^)0) + logZr is the dominant 
contribution, again because stabilizing selection affects the trait diversity less than its mean. 

5.3 Predictability of evolution 

How informative are trait measurements in one population about the distribution of trait 
values in a replicate population evolving in the same fitness landscape? To answer this 
question, we compare the ensemble-averaged Shannon entropy of the phenotype distribution 
within a population, 

(5)w = f S{W) QiW) (83) 
Jw 

and the Shannon entropy of the "mixed" distribution 

S{{W))=S{ f WQ{W)), (84) 
Jw 

which is obtained by compounding the trait values of all populations into a single distribution. 
We define the phenotypic predictability 

V = exp[{S)y^-S{{W))] (85) 

with S(W) = - J W{E) \ogW{E)dE. This quantity measures how much of the total trait 
value repertoire of all populations is already contained in the trait distribution yV{E) of a 
single distribution. It is closely related to the expected overlap between the distributions 
Wi{E) and W2{E) of two replicate populations. 

^Fitness flux plays a central role as a measure of adaptation also in non-equilibrium processes, where it 
is no longer related to free energy changes [41]. 
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To compute the predictability under stabilizing selection, we approximate the ensem- 
ble average in (83) and (84) by an average over F, using the approximate parametrization 
>V(^|r) - exp[-(E - r)V2(A)]. We obtain 

(A) / 1 



with the dimensionless ratio 

^ ^ {f)/{f)o ^1 [1 + 2c (1 + Q{ec))]~^ ^ (no recombination), ^^^^ 
(6)/ {S)o 1 [l + 2c (l + Qfree{c/i, «))] (free recombination) 

given by eqs. (64), (68), and (69). The dependence of V on the strength of stabilizing selection 
is shown in Fig. 5(c). While the neutral predictability Vo = 4^/(1 + 4:6) is small, stabilizing 
selection can generate predictability values V of order 1. The reason is again because the 
trait mean is more constrained than the trait diversity. This feature is illustrated in Fig. 1: 
under selection, a single-population distribution yV{E) fills a larger fraction of the trait range 
spanned by the cross-population distribution Q{T) than at neutrality. 

It is instructive to compare the phenotypic predictability (85) with the analogous measure 
for genotypes. 



Vg = exp[{S).,-Si{x) 

= exp [Xla ( (^a log Xa) - {x^) log{x^) ) ] 

For complex traits (i.e., for large values of £), we find the genotypic predictabihty 

P,~exp[-£[^(c)-0(^,ri)]]. (89) 
The leading entropy density <;{c) is given by 

^ ^ ^ - \ a E*/e -Jde K(e) log(l + e"^) for c > 1, ^ ^ 

where K,{e) is the single-locus effect distribution defined in eq. (70). The constant a is 
implicitly determined by the condition 



/ 



e e"^ E* 

de<e)-^ = -. (91) 



To derive this result, we note that q{c) is determined by the entropy of the "mixed" distri- 
bution, S{{x)), which can be evaluated in the low-mutation limit 9^0. Hence, ^(c) is also 
independent of recombination, which affects the overlap statistics between genotypes within 
a population [50] and only enters the ^-dependent corrections. Asymptotically for ^ <S 1 
and c <C 1, the mixed entropy reduces to the logarithm of the number of sequence states at 
the constitutive sites, S{{x)) ~ i log 2. In the strong-selection regime, we can compute this 



28 



entropy using the canonical formalism of statistical mechanics. We evaluate the partition 
function under linear selection on the trait, 

^ff = n 5Z = 11(1 + e'^^*) (92) 

i=l (Ti=0,l 1=1 

with the strength parameter a chosen to maintain the trait average at the fitness optimum, 

1=1 

The canonical entropy is then given by 5" = a{E) — log Zg — aE* — logZg, which leads to 
the entropy density (90). 

We conclude that the genotypic predictability is always small for complex traits, because 
an extensive number of genotypes remains compatible even with a strongly constrained 
trait value. Only after the projection from genotype to phenotype, selection can generate 
predictability. 



6 Inference of stabilizing selection 

Our results suggest a method to infer selection on a quantitative trait. The method is based 
on trait measurements within and across populations, but it does not require knowledge of 
the trait's genomic basis. Specifically, the ratio 

^,,,((r-(r)»^,,((r.-r.» 

(A) (A) 

depends only on phenotypic observables: it can be evaluated from the average trait diversity 
within populations, (A), and the variance of the trait mean across populations, ((F — (F))^), 
at evolutionary equilibrium (we assume the neutral sequence diversity 6 to be known inde- 
pendently). The ensemble variance ((F — (F))^) is just half of the equilibrium divergence, 
((Fi — F2)^), which, in turn, is close to the divergence between evolutionarily related popu- 
lations, ((r(ti) — r{t2)y), provided their divergence time is larger than the relaxation time 
of the trait to equilibrium. This is a reasonable approximation for traits under substan- 
tial selection, and our model can be extended to divergence data between closely related 
populations [28]. 

Our theory provides an analytical expression for Q, which is given by eq. (87). It shows 
that is a monotonically decreasing function of the strength of selection, c. This depen- 
dence can be used to infer c, which is defined as the fitness drop per 2N generations at 
a distance of one neutral standard deviation from the trait optimum. Both Q and c are 
pure numbers, which are independent of the units of trait and fitness. As shown by eq. (87), 
our phenotype-based method is formally similar to the well-known McDonald-Kreitman test. 



29 



0.1 05 1.0 5.0 10.0 50.0 100.0 

selection c 



Figure 6: Inference of stabilizing selection. The phenotypic observable Q measures the ratio 
between divergence and diversity of a quantitative trait, as given by eq. (94). Tliis ratio is plotted 
against the strength of stabilizing selection, c, for populations with different numbers {I) and effect 
distributions (k) of the trait's constituent sites, and with different recombination rates, (a) Data 
for non-recombing populations with i = 20, 100, 200 (dark to light green symbols) and two different 
effect distributions: delta distribution (all sites have equal effect, circles), exponential distribution 
(squares). Other system parameters as in Fig. 3. These data are in good agreement with the 
universal theoretical behavior Q{c) (solid line) given by eq. (87). Data points are shown within the 
range of applicability of the theory, c/i < 1 (for larger values of c, selection becomes balancing for 
individual loci), (b) Data for populations with free recombination for the same values of i (dark 
to light blue symbols) and the same effect distributions. These data are in good agreement with 
the theoretical behavior il(c) (lines) given by eq. (87), which contains a small dependence on i 
(dark to blue lines) and on the effect distribution (delta: solid lines, exponential: dashed lines), 
(c) Data for populations with different recombination rates p = 0.001, 0.01, 0.1, 0.5, oo (blue to green 
circles), evaluated for I = 100 and exponential effect distribution. These data interpolate between 
the theoretical predictions without recombination (green line) and with free recombination (blue 
line). Together, this shows the nearly universal dependence of the divergence-diversity ratio on the 
strength of stabilizing selection. 
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which evaluates divergence and diversity of genomic sequences [36]. However, the McDonald- 

Kreitman test has a different scope, which is to infer positive selection. 

The Q test exploits a universal characteristic of stabilizing selection: it affects the trait 
diversity less than its mean. This characteristic is quite intuitive from Fig. 1, which suggests 
that selection acts on divergence and on diversity with different characteristic strength. This 
strength is given by the curvature of the fitness landscape, Cq, multiplied with a relevant 
squared trait scale at neutrality. The basic such scale is the neutral expectation value of the 
trait divergence, (r2)o ~ E^. The trait scales within a population are different: Without 
recombination, selection acts on genotypes, and the relevant scale is the total trait diversity, 
9Eq. With strong recombination, selection acts on individual trait loci, and the relevant 
scale is the squared trait amphtude of one such locus, which is of order E^/i. Both within- 
population scales are small against the divergence scale E^. 

Most importantly, the inference of selection is confounded neither by number i and effect 
distribution k of the trait's constituent sites, nor by recombination between these sites. All 
of these genetic factors affect Q only through the term Q in eq. (87), which is small in the 
relevant range of 9 (at most percent) and i (at least tens of sites). As a result, Q depends on 
the strength of selection in a nearly universal way. Numerical simulations of populations with 
different site numbers, effect distributions, and recombination rates confirm this behavior, 
as shown in Fig. 6. 

7 Discussion 

In this paper, we have developed a statistical model for the evolution of complex molecular 

traits. We have shown that the dynamics of such traits can be described by approximate 
Kimura diffusion equations. In an arbitrary fitness landscape, this dynamics leads to cou- 
pled evolutionary equilibria for trait mean and diversity. Unlike the standard low-mutation 
or high-recombination approximations, our model is applicable to correlated multi-site pro- 
cesses, which evolve large hnkage disequihbria between the trait's constitutive sites. Such 
processes govern the evolution of complex traits in asexual populations; in sexual popula- 
tions, they are relevant for mesoscopic traits, which are polymorphic and based on a genomic 
region with limited recombination. Our model is a starting point for the analysis of such 
traits beyond the infinite-recombination assumption of quantitative genetics. It can and 
should be extended in a number of directions, which include the crossover from genotype 
selection to allele selection for finite recombination rates [50], traits with a nonlinear depen- 
dence on genotype, more rugged fitness landscapes, and time-dependent fitness "seascapes" 
driving adaptive trait evolution [40]. 

Our model leads to a new, quantitative test for stabilizing selection on quantitative traits, 
which is based on the ratio between trait divergence and trait diversity at equilibrium. We 
have shown that this ratio measures the strength of stabilizing selection in a nearly universal 
way, independently of the trait's genomic basis and of the recombination rate. This test can 
also be extended to quantitative traits in a time-dependent fitness seascape, which will be 
the subject of a forthcoming companion paper [28]. 
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Complex phenotypes integrate the information of multiple genomic sites. Compared to 
their constitutive genotypes, they represent biological functions on a larger scale. Both at 
the genomic and at the phenotypic level, we can ask about the predictability of evolution: 
How informative is sequencing or trait measurements in one population about the same 
quantities in a different population that evolves in the same fitness landscape? As we have 
shown in section 5.3, this question can be made precise by defining predictability in terms 
of an entropy difference between intra- and cross-population distributions of genotypes or 
trait values. For complex traits, predictability turns out to depend on scale and on selection. 
There is little predictability at the genome level, because the total number of genotypes 
encoding a functional trait is vastly larger than that realized in any one population. The 
equilibrium predictability is exponentially small in the number of trait sites, and populations 
evolving from a common ancestor will diverge through mutations at different sites. At the 
phenotypic level, the equilibrium predictability is related to the divergence-diversity ratio 
fl, as given by eq. (86). It is small at neutrality, but under sufficiently strong stabilizing 
selection, it can reach values of order one. Hence, stabilizing selection generates predictability 
of evolution at the phenotypic level. 
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Appendix: 

Non-equilibrium ensembles of quantitative traits 

Here we analyze the evolution equation (35) for the joint distribution Q(r,A,t) of trait 
mean and variance in asexual populations. We consider the case of stabilizing selection in the 
fitness landscape F{T, A) given by eq. (58). Using the scaled trait variables 7 = {T — To)/Eo, 
5 — A/Eq and A = (F — E*)/Eq and the scaled selection strength c = 2NEqCo defined in 
eqs. (59) and (61), this equation can be written in the form 



d_ 

d_ 
dS 



^6 + 46^ + 2c 5(7 - (A)o) 
d 



dS 



2S^ + 89{S - 1) + 25 + 2c 



(95) 



where {J'^, J^) denotes the probability current in the ^-S plane. The scaled distribution and 
current are related to their unsealed counterparts, 



Q{T,A-co,Eo) 
J^(r,A;co,^o) 
J^(r,A;co,£;o) 



E-'/'Q{^,6;c) 
E^'r{^,S;c), 
E-'/'j\^,d;c) 



(96) 
(97) 



If we keep the trait effect distribution «;(e) of individual constituent sites fixed, the squared 
overall trait scale is proportional to their number, E^^ = K2i'/4. Hence, we can interpret the 
relations (96), (97) as scale transformations relating systems with different values of £. 

The joint dynamics of F and A does not have an equilibrium solution, because the muta- 
tion coefficient field (4^7, S9{6—l)+26) is non-integrable. As shown by numerical simulations, 
this dynamics leads instead to a non-equilibrium stationary distribution Qstat(7;^) (shown 
in Fig. 2), which has with a finite current {J'^ , J^){^,6). According to eqs. (96) and (97), 
we obtain a scale-invariant stationary non-equilibrium state^ describing a family of systems 
with different £ and constant scaled selection strength c. In this family of systems, a finite 
stationary current persists for large values of £. 

To test this prediction, we measure the current by binning changes of F and A through 
discretized grid lines in our simulation. At each junction in the grid, we record positive and 
negative changes separately. These changes determine the current components J'"(F, A) = 
jr(r,A) - J^(F,A) and J'^(r,A) = J^(r,A) - (r,A), and we obtain dimensionless 
measures for the violation of detailed balance, 

J^(F) - J^(F) 



Jf(F) + J^(F)' 



J^(A)-J^(A) 
jr(A) + jr(A)' 



(98) 



^The scale-invariant state (96, 97) of asexual neutral evolution (c = 0) can be broken by two relevant 
perturbations: recombination (which changes the scaling of diversity fluctuations from A ~ £ to A ^ fi^^) 
and stabilizing selection of constant unsealed strength cq (which generates a trait autocorrelation (C) < 0). 
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Figure 7: Deviations from equilibrium under complete genetic linkage, (a) Current ratios 
Jsij) (left panel) and J'y{6) (right panel), as defined in eq. (98). (b) Deviations R{'^) = EqR{T) (left 
panel) and R{6) = EqR{A) (right panel) of the stationary distribution from the Boltzmann form, as 
defined in eq. (99). These data are shown for systems with i = 100, 120, . . . ,200 constitutive sites 
of equal effect. They collapse into unique functions, indicating a scale-invariant non-equilibrium 
state with stationary current; see eqs. (96) and (97). The fitness optimum is set to E* = Fq = i/2. 
Other system parameters are as in Fig. 3. 

which are obtained from the integrated currents t7±(r) = J J^(T, A) dA and J±{A) = 
J ^±(r, A) dr (these measures are less noisy than their local counterparts). In Fig. 7(a), we 
show the ratios (98) for systems with different values of £ (at a fixed value of c). As predicted 
by the scale transformation (97), these data collapse onto unique functions J'si'y) and J^y{6). 
They reveal a closed stationary current with a clockwise loop for 7 > and a symmetrical, 
counterclockwise loop for 7 < 0. 

The breakdown of detailed balance in the stationary state has an important consequence: 
the distribution (5stat(r) A) under selection is no longer of Boltzmann form. The deviations 
are measured by the function 

R{T, A) = Q,tat(r, A) - ^ Qo{T, A) exp[2iVF(r, A)]. (99) 

Fig. 7(b) shows the marginal differences R{T) = J R{T, A)dA and R{A) = J i?(F, A)dT for 
systems with different values of £ (at a fixed value of c). After rescaling according to eq. (96), 
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these data again collapse onto a single function R{^,S). The actual stationary distribution 
<5stat(r, A) is seen to be broader than the corresponding Boltzmann distribution. 
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